Master equation approach to computing RVB bond amplitudes 
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We describe a "master equation" analysis for the bond amplitudes h{r) of an RVB wavefunction. 
Starting from any initial guess, ^(r) evolves — in a manner dictated by the spin hamiltonian under 
consideration — toward a steady-state distribution representing an approximation to the true ground 
state. Unknown transition coefficients in the master equation are treated as variational parameters. 
We illustrate the method by applying it to the J1-J2 antiferromagnetic Heisenberg model. Without 
frustration (J2 = 0), the amplitudes are radially symmetric and fall off as 1/r^ in the bond length. 
As the frustration increases, there are precursor signs of columnar or plaquette VBS order: the bonds 
preferentially align along the axes of the square lattice and weight accrues in the nearest-neighbour 
bond amplitudes. The Marshall sign rule holds over a large range of couplings, J2/J1 < 0.418. It 
fails when the r — (2, 1) bond amplitude first goes negative, a point also marked by a cusp in the 
ground state energy. A nonrigourous extrapolation of the staggered magnetic moment (through this 
point of nonanalyticity) shows it vanishing continuously at a critical value J2/J1 ~ 0.447. This may 
be preempted by a first-order transition to a state of broken translational symmetry. 



I. INTRODUCTION 

In the early 1970s, a resonating-valence-bond (RVB) 
wavefunctioni with nearest-neighbour (NN) bonds only 
was proposed as a possible ground state for the quantum 
Heisenberg model on the triangular lattice^i^ This short- 
ranged, quantum-disordered^ RVB state was conceived 
in analogy with the spin liquid state found in one di- 
mension.^'^ The belief was that classical 120° Neel order 
was unlikely to survive in the presence of strong quantum 
fluctuations. 

This conjecture ultimately proved incorrect. Like other 
low-coordination-number antiferromagnetSf^ the trian- 
gular system is ordered at zero temperature J^i^iii^ Con- 
sequently, its ground state cannot be described in a ba- 
sis of short bonds. One can show, in fact, that a cor- 
rect description must involve valence bonds on all length 
scales 

A generalization of the RVB state that includes long 
bonds was later proposed by Liang, Doucot, and Ander- 
son for use as a variational wavefunction in the square- 
lattice Heisenberg modelJ^ Their idea was to factor- 
ize the weight associated with each valence bond con- 
figuration into a product of individual bond amplitudes 
that depend only on the vector r connecting bond end- 
points. Unlike the NN-bond RVB, which is unique, the 
long-range version is a family of states parameterized 
by the bond distribution function, h{r). In d = 2, the 
RVB wavefunction has expressive power to describe both 
an antiferromagnetically ordered phase and a featureless 
quantum disordered phasc;^-'^^ It may be a good vari- 
ational wavefunction for systems in which antifcrromag- 
netism is killed by the addition of frustrating interactions. 

As a practical matter, optimizing the bond amplitudes 
numerically is not straightforward. The number of inde- 
pendent parameters is of the order of the system size, and 
the energy depends only very weakly on the amplitudes 
of the longest bonds. Thus, obtaining well-converged re- 
sults becomes increasingly difficult for large lattices, and 



scaling to the thermodynamic limit is unreliable. Lou 
and Sandvik-lS have made some progress by experiment- 
ing with different optimization schemes. They recently 
carried out an unbiased variational determination of h{r) 
for the square-lattice Heisenberg model and were able to 
achieve lattice sizes up to 32 x 32. 

Liang, Doucot, and Anderson circumvent the problems 
associated with a macroscopic number of degrees of free- 
dom by assuming a functional form for h(r). They vary 
the amplitudes of only a few short bonds and fix the re- 
mainder under the assumption of a radially symmetric 
bond-length distribution and algebraic decay at long dis- 
tances,^^ For local, nonfrustrating interactions, this as- 
sumption turns out to be essentially correct.^'' Nonethe- 
less, their choice of functional form is ad hoc, and there 
is nothing in their approach that provides insight into 
how the amplitudes might change when competing inter- 
actions are introduced. 

In this paper, we describe an alternative method for 
calculating the bond amplitudes that requires at most a 
few variational parameters. The utility of the method is 
tested by applying it to the J1-J2 model. As in Ref. [3, 
we make strong assumptions about the form of the bond 
distribution. In our case, however, the choice of func- 
tional form for h{r) is guided by a master equation that 
mimics the reconfiguration of bond amplitudes induced 
by the evolution operator. 

The J1-J2 model describes a system of spin-half mo- 
ments arranged on a square lattice in which Heisenberg 
interactions of strength Ji, acting along the plaquette 
edges, compete with frustrating interactions of strengh 
J2 , acting across the plaquette diagonals. At J2/ Ji = 
and J2/J1 = 00, the model has two- and four-sublattice 
Neel order, respectively. There is a gapped intermediate 
phase in the vicinity of J2/J1 ~ 0.5, whose exact nature 
remains controversial. There has been speculation about 
a possible spin liquid statep2ii2S'2iiS2i22- but a state with 
broken translational symmetry now seems more likely. 
The leading candidate is a valence bond solid (VBS) with 



either columnar— or plaquett o^^'^^i'^° order. 

The extent of the intermediate phase has been deter- 
mined to about one digit of precision. Exact diagonal- 
ization on small clusters^ puts the lower critical point 
at J2/ Ji = 0.34(4), but this appears to be an under- 
estimate. Bond operator calculations^Si^ based on the 
columnar VBS predict 0.38 < Jij J\ ^ 0.62 for the region 
of stability, and series expansion s'^ "^'^^ from the magnetic 
side give 0.4 < J2/J1 < 0.6. A quantum Monte Carlo 
studyr^ in which stochastic reconfiguration is used to 
partially alleviate the sign problem, reports a transition 
to a gapped state at J2/J1 » 0.4. 

It has been established from energy level crossings in 
series expansion that the transition at the upper criti- 
cal point is first order^^^^ No such crossings have been 
detected at the lower critical point, at least within the 
numerical accuracy that can be achieved. In most of 
the studies cited above, it is implicitly assumed that 
the transition at the lower critical point is second order. 
If that is true — and if the intermediate phase is indeed 
bond ordered — then the lower critical point may consti- 
tute a deconfined quantum critical point, as envisioned 
by Senthil ei ali^ 

The fact that bond operator methods indicate a high 
density of triplet modes near a deconfined quantum criti- 
cal point,— but only a low density near the critical point 
of the J1-J2 modelr^ leaves room for doubt. Indeed, a 
recent series expansion study points to a first order tran- 
sition at J2/J1 ~ 0.43 on the basis of an energy func- 
tional computed for a fictitious translational-symmetry- 
breaking field. This is supported by the argument due 
to Chubukov^^ that a continuous transition is only pos- 
sible when a third-nearest-neighbour interaction J3 > 
is present. 

The results reported here cannot settle this question 
with any certainty, but they do appear to be more con- 
sistent with a first order Neel-VBS transition. 



II. MASTER EQUATION FOR FACTORIZABLE 
RVB BOND AMPLITUDES 

The spin-rotation-invariant (total spin S = Q) ground 
state of a system of 2N spin-i moments can be written 
as a superposition of valence bond states42ii^ The sim- 
plest RVB ansatz is to assume that the weight given to 
each bond configuration is a product of individual bond 
amplitudes: 




FIG. 1: Spins on the A sublattice (solid circles) and B sub- 
lattice (open circles) are grouped into pairs forming singlets. 
Each pairing configuration is characterized by a set of directed 
bonds connecting A sites to B sites. 



whenever r connects valence bonds in the same sublat- 
tice. See Fig. [TJ This is just a restatement of the Marshall 
sign theorem."*^ 

One way to compute the h{r) values appropriate for 
a given model is to consider the r-dependent family of 
states 



\Kr)) = 



IMO)), 



(2) 



where H is the hamiltonian of interest and T is an op- 
erator that projects onto the space of factorizable RVB 
wavefunctions. In each time step dr, some fraction of 
the bond amplitude is reapportioned as bonds are cre- 
ated and destroyed. Correlations between bonds that go 
beyond the RVB framework are prevented from accumu- 
lating. This process is governed by a master equation 
that describes how the distribution h{v) evolves towards 
its steady-state solution. Note that the wavefunction 
that emerges in the r — > c» limit is not strictly equal 
to the projection T\i>) of the true ground state 1-0); nor 
is it equal to the variationally determined state that 
minimizes E — {h\H\h) / {h\h) . Nonetheless, all three are 
very similar to one another4^ 

The key observation is that the valence bond basis 
is closed under operation by the Heisenberg interaction. 
Operating on an existing bond simply leaves the bond as 
is [and the distribution h{v) unchanged] whereas operat- 
ing between two bonds maps them to their complemen- 
tary tiling: 



(3) 



(1) 



i-S,.S,)[z,Z][fc,j] = i[z,j][fc,Z]. (4) 



Here [i, j] = ^(||j|j) - |ijtj) denotes a singlet formed 
from the spins at sites i and j. The effect of Eq. ^ is 
depicted in Fig. [51 

For the NN Heisenberg model on a d-dimensional 
(hyper-)cubic lattice, the master equation is 



Here the sum is over all partitions of the lattice into N 
singlet pairs, and the product is over all vectors r drawn 
between valence bond endpoints. [Anderson's NN-bond 
RVB corresponds to /i(r) = 5[\r\ — 1).] In the special case 
of a nonfrustrated model on a bipartite lattice, the am- 
plitudes h(v) are real and nonnegative and strictly zero 



h{v) = ^K,a+5I <5r'+r"-a,r/i(r')/i(r")]-2z/i(r), (5) 



where h = dh/dr, z — 2d is the coordination number, 
and a ranges over all NN vectors. This is correct only 
insofar as h(r) accurately measures how often a bond of 
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FIG. 2: A nonfrustrating Heisenberg interaction, indicated 
by the dotted (red) line, is applied between sites in opposite 
sublattices. The resulting reconfiguration creates one valence 
bond where the interaction was applied and another between 
the two remaining endpoints. 



type r appears in the superposition of valence bond con- 
figurations making up the RVB state. Geometrical tiling 
constraints, which are increasingly important at low co- 
ordination number, have been ignored. Nonetheless, this 
level of approximation allows us to proceed analytically. 

Equation ([5]) conserves the unit normalization of the 
total weight: 
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FIG. 3: The bond amplitude functions predicted for the 
linear-chain (upper, blue curves) and square-lattice (lower, 
green curves) Heisenberg models are plotted for various sys- 
tem sizes. Larger values of L are indicated by darker data 
points, following the legend. The solid lines reflect the ana- 
lytical results given in Eqs. (|10|l and (lll|l . 



z + z 



J2hir) -2zJ2hiv)^0. (6) 



The AB character of the bonds is also a constant of 
the motion. If we start with a distribution h{r) that 
is nonzero only when r connects sites in opposite sublat- 
tices, then h{r) will also have this property at all subse- 
quent r. 

A somewhat stronger property of the flow is that 
all weights associated with bonds of even Manhattan 
length ||r|| = |ri| -|- \r2\ + ■■■\rd\, namely the AA or 
BE bonds, are driven to zero. This is a straightfor- 
ward consequence of an asymmetry in the reconfiguration 
rules: (even, odd) (even, odd), (odd, odd) — > (odd, odd), 
and (even, even) —!■ (odd, odd). This is yet another mani- 
festation of the Marshall sign rule. 

Accordingly, for r — > oo there are no bonds connecting 
sites in the same sublattice and all bonds have odd Man- 
hattan length. We are thus free to impose the convention 
that the vector character of all bonds is directed from A 
to B (as anticipated in Fig.[T]). This means that the bond 
amplitude function has a Fourier expansion 



(7) 



where the wavevector sum ranges over a reduced "mag- 
netic" Brillouin zone, equal to the standard Wigner-Seitz 
cell modulo Q = (tt, . . . tt). One finds that the Fourier 
transform of Eq. ([5]) is a simple polynomial in h^, 



q' 



whose stationary distribution is 



/iq = 



2)1/2 



7q 



(8) 



(9) 



7q = (l/d)(cosgi cos qd) is the Fourier transform of 
the NN matrix. In real space, the long distance behaviour 
is given by 



h{r) 
h{r) 



7r(l + r2) 
V2 

27r(i -f-r2)3/2 



(d=l) (10) 
(d = 2) (11) 



as shown in Fig. [31 Note that in two dimensions, the bond 
amplitude is almost prefectly radially symmetric beyond 
a few lattice spacings. The general behaviour for higher 
dimensions is h{r) ~ {1/d + r^)-'^'^+^^/^ ~ r'^d+il 



III. FRUSTRATING INTERACTIONS 

As we emphasized in the previous section, any model 
on a bipartite lattice whose interactions are nonfrus- 
trating with respect to two-sublattice Neel order can 
be described in a basis consisting only of AB valence 
bonds.— i^i^ Two special features of the AB basis are 
that (1) the overlap between any two states is strictly 
positive;^ and (2) there is an exact correspondence be- 
tween the Marshall sign rule and the positivity of all the 
RVB bond amplitudes. 

We now argue that, even with the addition of frustrat- 
ing interactions, one can still chose to work exclusively in 
the AB basis. According to Eq. ([3]), a frustrating inter- 
action applied between sites in the same sublattice trans- 
forms two AB bonds into one AA and one BB bond. But 
since valence bonds are nonorthogonal, we can take ad- 
vantage of the overcompleteness relation 



(12) 



4 




FIG. 4: A frustrating Heisenberg interaction applied between 
two A sublattice sites, indicated by the dotted (red) line, has 
the effect of exchanging the two valence bond endpoints. 



to eliminate each of the unwanted bonds, yielding a new 
update rule 

(l + S,.Sfc)[z,;][A:,j] = i[z,j][fc,;], (13) 

where i,k £ A and j,l S B. See Fig. [H There is no 
diagonal operation analogous to Eq. 

Following Eq. , a model with NN Heisenberg inter- 
actions of strengh Ji and next-nearest-neighbour (NNN) 
interactions of strength J2 has the bond amplitude mas- 
ter equation 



r 4-r 



,h{r')h{r")] - 2zh{r) ■ 



J2 
Ji 



'Jr"+a.,r 



h{r')h{r")-2Sh{r)], (14) 



a r ,r 



where a ranges over all NNN vectors. (We use a tilde to 
distinguish NNN quantities from NN ones.) This differs 
from Eq. ([5]) by a term proportional to J2/ Ji- 
Fourier transformation of Eq. leads to 



/iq = 7q + 7q^q - 2 



l+5(l-7q) 



hr 



This has a steady state solution 



A, 



7q 



where 



Aq = 1 + 5(1 - 7q) and g ^ {z / z){ J2/ Ji). 



(15) 



(16) 



(17) 



In dimension d > 1, the coordination numbers are z = 
2d, z — 2'^ and the connection matrices have Fourier 

transforms 7q = ^ Ylk=i (^osqk and 7q = Jlfc^i cosgfe. 

Figure [5] illustrates the real-space distribution corre- 
sponding to Eq. (|16p for several values of g in two di- 
mensions. When g < 0, the long-range behaviour is 



hir) 



\ as in Eq. ([TT|) . When g > 0, the radial 



symmetry is reduced to the C4 symmetry of the square 
lattice, and the amplitudes begin to accumulate along the 
principle axes, especially in the r — (1, 0) bond. Contrary 
to our expectations, the distribution does not become 
uniformly more short-ranged as g increases. Along the 
principle axes, it actually becomes longer-ranged: the ex- 
ponent of the algebraic decay steadily decreases from 3 
(at 5 = 0) to 1.5 (at g = 0.5). This looks nothing like 
the h{r) ~ r^P spin liquid found at p > 3.3.^"''^^ 

The Marshall sign rule is obeyed in the J1-J2 model 
up to relatively large values of the frustration param- 
eter4^>^>^ At the level of approximation employed 
here, the bond amplitudes are all positive up to gM = 
0.323158, the coupling at which the amplitude of the 



r — (2, 1) bond passes through zero. A sign change in 
h{2, 1) at large frustration has also been observed by Lou 
and Sandvik in their unbiased calculation.^^ 

Since there is already some ambiguity in the master 
equation because of the neglect of geometric constraints, 
we will treat g as a variational parameter. In other 
words, we will allow the relative weighting between the 
frustrating and nonfrustrating channels to deviate from 
g = (z/z)(J2/Ji), as the energy dictates. 

In principle, the variational approach can be expanded 
to include farther- neighbour moments, defined by 



7q^ = 



1 

^ E n cos(^<,(„)fc„). 



(18) 



<y£Sd n=l 



Here, the index fi is an ordered d-tuple of natural num- 
bers, and Sd is the set of permutations on d elements. 
There will be a variational parameter g'^ for each in- 
cluded moment, in terms of which the amplitude distri- 
bution is 



[(A-ryq)2-A2] 



1/2 



Ac 



where 



^q = E«' 



and 



Aq = E«- 

odd 



(19) 



(20) 



A is fixed by the requirement that ft.q=o = 1- The sum- 



mations in Eq. Ij20p are over all fx vectors having even 
and odd Manhattan length up to some cutoff, ||/x|| < fiQ. 
For our numerical work on the J1—J2 model, only the 
fi = (1,0) and fjb = (1, 1) components are kept. 
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FIG. 5: The bond amplitudes h{r) are depicted from left to right for the values g ^ 0, 0.2, 0.32, 0.45, 0.49999. The top row 
shows the subset of bonds up to length (16,15) for an L = 256 square lattice. The bottom row is a magnified view emphasizing 
the short range bonds up to (6,5). The position of each circle marks the endpoint of a bond whose other endpoint is at the 
origin. The area of each circle is proportional to \h{r)\r^^^ . Filled (black) circles denote a positive value and open (red) circles a 
negative one. For <? < 0, the bond amplitudes are radially symmetric and positive definite and fall off as r"''. As g increases, the 
distribution becomes increasingly asymmetric. The (2,1) bond steadily decreases in magnitude and vanishes at qm = 0.323158. 
In the limit g 0.5, the bonds along the x and y axes become extremely long ranged: h{r,Q) = h(0,r) ~ r~^^^. For g > 0.5, 
the bond amplitudes are complex. 



IV. RESULTS FOR THE J1-J2 MODEL 

We work with an RVB trial wavcfunction whose 
weights are factorizcd as in Eq. The bond ampli- 

tudes are taken from the Fourier transform of Eq. ()16|1 . 
These depend only on the size of the lattice and on a 
single variational parameter, g, which is fixed by mini- 
mizing E{g) = (H). Expectation values of an operator 
O in the trial state, written 



{v\v') 



(21) 



can be interpreted as an ensemble average of the estima- 
tor {v\0\v') / {v\v') in a fluctuating gas of valence bond 
loopsi^i^ Over the range g < gM, the bond amplitudes 
are all strictly positive and thus the sampling weight 



{v\v 



l[h{r) llh{r' 



(22) 



has no sign problem associated with it. Numerical eval- 
uation of the RVB wavcfunction is carried out using 

a worm algorithm^^ adapted to the valence bond loop 
54 

gasi^ 

FigureElshows the NN and NNN spin correlations com- 
puted for finite lattices as a function of g. These data 
are extrapolated to the thermodynamic limit by assum- 
ing 0{L^^) leading corrections. A weighted sum of the 
correlations gives the variational energy: 



Eig) 



2^ 



E 



-(1,0); 



J2 
Ji 



r+(l,l)/ 



(23) 



The optimal value of g is found by solving E'(g^i^) = 0. 
The dependence of gmin on J2/J1 is shown as an inset 
in the bottom panel of Fig. [6l Back substitution of ^min 
into Eq. 1^ gives E{J2/Ji). 

As is clear from the top and middle panels of Fig. [51 
a point of nonanalyticity at g — gM (or J2/J1 — 0.418) 
separates regions with markedly different behaviour. In 
the L 00 limit, both the NN and NNN spin correla- 
tions exhibit a cusp. E'{g) has no roots for g > gM- 

Figure [7] shows the staggered magnetic moment 



M 



1 

2N 



-,1/2 



(Sr • Sr') 



(24) 



plotted as both M{g) and Af(J2/Ji). Here, the ther- 
modynamic limit is acheived with 0{L~^) scaling. By 
continuing the trend established in the region where the 
Marshall sign rule is obeyed, we estimate that the stag- 
gered moment vanishes continuously at a critical coupling 
J2/J1 = 0.447. 



V. DISCUSSION 

The master equation approach applied to the J1-J2 
model leads to an RVB trial wavcfunction whose bond 
amplitudes depend on a single variational parameter, g. 
Wherever the Marshall sign rule holds, we are able to 
compute the properties of the RVB state to very high 
accuracy for large lattices (up to L = 128 easily on 
a laptop) and thus to extrapolate measured values to 
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FIG. 6: (Top) The expectation value of the NN spin correla- 
tions iSi = —{1/2N) X]r(^r ■ Sr-i-(i,o)) IS plotted as a function 
of the variational parameter g. Solid circles (blue) represent 
data computed for a particular L x L system; see the legend. 
Error bars (red) denote the L — > 00 extrapolation. (Middle) 
The NNN spin correlations ^2 = (l/2iV)Er(Sr • Sr+(i,i)) 
are plotted in the same way. (Bottom) The energy density 
E/ Ji = —Si + {J2/ J\)S2 in the thermodynamic limit is com- 
pared to estimates (black dots) due to Gochev.^^ The inset 
shows the g that minimizes the variational energy. 



0.5 



0.3 



M 



1 1 


— 1 1 1 1 


— 1 1 1 


^^^^^wiiiiiiiiiiiiii 


'^i^iii;!:;:;:!!^:;;::::''" 










• 8 


• 32 




• 12 


• 48 




• 16 


• 64 




- • 20 


• 96 




• 24 

■ I 


• 128 
■ III 





-0.4 



-0.2 



0.2 



0.4 



M 



0.5 
0.4 
0.3 
0.2 
0.1 




-0. 




• 8 


• 32 


• 12 


• 48 


• 16 


• 64 


• 20 


• 96 


• 24 


• 128 





0.1 



0.2 
J2/J1 



0.5 



FIG. 7: (Top) The staggered magnetic moment M is plotted 
as a function of the variational parameter g. The dots (in 
shades of blue) represent data computed for a particular LxL 
system. Errorbars (red) indicate the extrapolated L ^ 00 
value. The data for g > gM is unreliable for reasons elucidated 
in the main text. The solid line is a polynomial fit of the 
L — 00, g < gM points. (Bottom) The g < gM data from 
the top panel is replotted with a new horizontal scale. The 
Marshall sign rule breaks down at 0.418, and the magnetic 
order dies at 0.447. 



the thermodynamic hmit. For the NN Heisenberg model 
( J2 = 0), the energy and staggered magnetization of the 
best variational state (at g = —0.0484) extrapolate to 
E = -0.669748 and M = 0.3086. These differ by 0.047% 
and 0.52% from the exact results E = -0.669437(5) and 
M — 0.3070(3) obtained from quantum Monte Carloi^ 

As a function of frustration, the antiferromagnetic or- 
der dies out quite slowly. When the Marshall sign rule 
finally breaks down (at J2/J1 = 0.418), the staggered 
moment has decreased only to M ~ 0.1114, about 36% 
of its unfrustrated value. This does not appear to be con- 
sistent with a continuous transition. One point of concen- 
sus for this model is that a gapped state appears around 
J2/J1 w 0.4. Thus, the fact that the magnetization is 
still large near this value suggests that the transition is 
first order. In comparison to other situations where a 
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Neel-VBS transition is known to occur what we ob- 
serve here is much more Uke the situation in Ref. [sgI . 
where the staggered magnetization decreases only mod- 
estly and then collapses abruptly at a first-order critical 
point. 

Of course, this line of reasoning is not sufficiently 
rigourous to establish the order of the transition, and 
we cannot rule out a deconfined quantum critical point. 
Since there are no exact results on large lattices, we do 
not know how well the optimized RVB wavefunction ap- 
proximates the true ground state. (For 4x4 the agree- 
ment is quite good^-: the overlap is 0.9998 for the un- 
frustrated case and 0.996 for J2/J1 = 0.4). We have esti- 
mated that the magnetic order vanishes at J2/ Ji — 0.447 
in a continuous scenario, which apparently "overshoots" 
the openning of the spin gap at J2I J\ ~ 0.4. It is hard to 
say whether these values are truly noncoincident, since 
their uncertainties are difficult to quantify. Moreover, the 
decay of the staggered magnetization may be artificially 
slow because of the failure of the RVB state (which is 
translationally invariant) to capture the incipient dimer 
correlations near the transition. 

The sign problem in this model (for J2/J1 > 0.418) 
turns out not to be terribly severe. Much more catas- 
trophic is that the master equation itself breaks down 
along with the Marshall sign rule, because of the assump- 
tion that /i(r) > represents the probability of finding 
a bond of type r. Once any of the amplitudes becomes 
negative, the reasoning that lead to Eq. is no longer 
correct. The breakdown could perhaps be avoided if we 
were to use an exact numerical implementation of Eq. ^ 
to find the r — s- cx) limit, rather than an analytical ansatz. 
More likely, though, the failure of the master equation is 
related to the inability of the RVB state to accommodate 
bond-bond correlations — except indirectly by strength- 



ening the C4 symmetry of /i(r), as seen in Fig. [5l 

The master equation approach works remarkably well 
in guiding our choice of the RVB bond amplitudes. 
Where it can be checked ( J2 = 0), the accuracy of the 
wavefunction rivals that of unbiased optimizations, but 
with an enormous computational saving associated with 
reducing the number of variational parameters from TV to 
1. Including variational parameters for a few additional 
modes (as described at the end of Sect. IIII[) would im- 
prove the accuracy further. In order to handle the most 
disruptive effects of frustrating interactions, however, it 
will be necessary to move to the next level of approxima- 
tion and to consider RVB states whose weights factorize 
into amplitudes for 'pairs of bonds. Obtaining an analyt- 
ical master equation for the two-bond amplitude hij-ki, 
as we did in this paper for the single-bond amplitude 
hij = h{rij), is probably not feasible. Nonetheless, for 
variational calculations, it may be enough to put in by 
hand some bond-bond contribution, e.g., 

hiy.ki = h[vij)h{vki) 

X [l + U5{\\v,,+Vki\m\\v,i+vu,\\)\, (25) 

that is compatible with the expected VBS pattern (here, 
the columnar state at large U). 
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